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Abstract 


A  method  based  on  an  augmented  Lagrangian  formulation  is  developed 
which  allows  one  to  estimate  coefficients  in  an  elliptic  differential  equation 
from  measurements  of  the  state.  This  is  a  hybrid  method  combining  the 
output-least-squares  and  the  equation-error  technique.  Seminorm  regular¬ 
ization  is  employed,  and  convergence  and  stability  properties  are  discussed. 
Several  aspects  of  an  efficient  implementation  are  decribed.  Finally  the  effec¬ 
tiveness  of  the  method  is  demonstrated  by  means  of  one  and  two  dimensional 
examples. 
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1  Introduction 


In  this  paper  we  discuss  our  numerical  experience  with  the  hybrid  method 
[IK2]  based  on  the  augmented  Lagrangian  formulation  for  estimating  the 
coefficient  q  in  the  elliptic  equation, 

(1.1)  — div  (g  gradu)  =  /  in  fi 

t/|r  =  0 

from  a  measurement  2  of  the  solution  u  (9’)  which  corresponds  to  the  “true" 
coefficient  9*.  Here  0  is  a  bounded  open  set  in  =  1,2,3  with  piecewise 
smooth  boundary  F,  and  9  is  a  scalar  valued  function  of  the  spatial  variable 
X.  Frequently  n{q’)  is  not  known  exactly;  this  may  be  due  to  model  and/or 
measurement  error  and  also  because  measurements  may  be  known  only  at 
discrete  points  i,  in  the  domain  Q  so  that  2  is  constructed  by  interpola¬ 
tion  of  pointwise  data.  A  motivation  for  developing  the  hybrid  method  is 
to  combine  the  output  least  squares  method  and  the  equation  error  method 
[C]  into  one  algorithm  while  retaining  the  favorable  properties  of  both.  The 
augmented  Lagrangian  formulation,  just  as  the  equation  error  method,  re¬ 
duces  to  a  quadratic  programing  problem  and,  like  the  output  least  squares 
method,  it  is  versatile  with  regards  to  availability  of  observations  as  will 
be  demonstrated  in  Section  4.  The  method  is  based  upon  formulating  the 
problem  as  a  constrained  and  regularized  minimization  problem  (P^)  for  two 
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independent  variables  (q,  u)  €  x  and,  as  will  be  shown  shortly. 

the  minimization  (P^)  is  solved  interatively  by  the  augmented  Lagrangian 

method.  Let  us  consider: 

(P^)  minimize  \\U-Z  i;,,  +fAr(,) 

e{q,  u)  =  (— A)“*  (drv  {qgrad  u)  +  /)  =  0 

subject  to  l{q)  =  a  —  q{x)  <  0 

5(9)  =  $(!  9  !//2  -7^)  <  0, 

where  we  assume  that  /  €  H~^{Q)  and  z  €  (H),  that  q  and  7,  are  positive 

constants  satisfying  o  meas(17)'^^  <  7  and  where  A  denotes  the  Laplacian 
with  Dirichlet  boundary  condition  as  an  operator  from  Hq  (fl)  onto 
Note  that  c  :  x  Hq{Q)  — ♦  Hq{Q)  is  a  continuous  bilinear  map.  and 

that  g  :  — *  3?’  and  the  affine  map  C  :  are  continuous. 

Finally  ^N{q)  is  a  regularization  term. 

It  is  easy  to  show  [IK2]  that  there  exists  a  Lagrange  multiplier  A'  €  Hq{9.) 
associated  with  the  equality  constraint  e{q,u)  =  0  and  it  is  given  by 

(1.2)  A- =  (v(/v))-' A 

where  {q^^.u^)  denotes  a  solution  of  (P"^). 

A  number  of  remarks  are  in  order  (  see  also  Section  3  for  a  further  dis¬ 
cussion). 

(1)  The  pointwise  constraint  £  guarantees  the  strong  ellipticity  and  the 
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norm  constraint  g  allows  us  to  argue  the  existence  of  solutions  of  (P^),  > 

0.  Without  the  norm  constraint,  (P^)  may  not  have  a  solution  in  general. 

(2)  The  second  term  in  the  objective  functional  represents  a  regularization 

term  in  which  >  0  and  N{0)  is  a  seminorm  on  e.g. 

(1-3)  A'(«)  =  1  li.  +  E  I Ii>  • 

l<»j<n 

For  an  appropriate  choice  o{  0  >  0  the  use  of  a  regularization  term  in  (P^) 
induces  the  continuity  of  the  mapping  from  the  observation  z  £  \o  the 
solution  (^(r),  ti'’(s))  of  (P^),  see  [CKl,  CK2].  In  general  ^  cannot  be 
taken  equal  to  zero  if  continuous  dependence  of  the  solutions  of  (P'^)  on  z  is 
I  desired. 

(3)  bince  u  is  uniquely  determined  from  q  through  the  equality  c(q.  it)  = 
0.  the  objective  functional  can  be  considered  to  be  dependent  on  q  only.  In 
this  way  (P*^)  is  reduced  to  the  well  known  regularized  least  squares  method 
given  by 

minimize 

subject  to  —  div(^  grad  u{q))  =  /,  l{q)  <  0.  and  g{q)  <  0. 

From  an  optimization  point  of  view  {P^)  and  the  least  squares  method  are 
equivalent.  But  in  (P^)  q  and  u  are  independent  variables  and  hence  the  im¬ 
plementations  based  on  (P^)  are  very  different  from  those  of  the  regularized 
least  squares  formulation. 


5 


(4)  If  rj  =  1  the  parameter  space  can  be  replaced  by  //’(fj).  In 

this  case  g{q)  =  |(l9l«i  -  h'^)  and  N(q)  = 

The  augmented  Lagrangian  method  applied  to  {P^)  involves  a  sequence 
of  minimizations  of  the  functionals 

<  .  i  I  u-z  j^,  +  §N(q) 

+  (^'‘,e(q,u)}ffi  +  ^  1  e(9.t/) 

subject  to  {(q)  <  0  and  ^(9)  <  0 

and  the  multiplier  sequence  {A^}  in  generated  by  the  updating  rule 

(1.5)  A‘+*=A‘  +  c,c(9^,Ufc) 

where  the  pair  (q^,  u*)  is  a  (local)  minimizer  of  lck(’^  s  .^i).  To  carry  out  this 
iterative  scheme  a  sequence  of  monotonically  nondecreasing,  positive,  real 
numbers  {cj,-}  and  the  start-up  value  XUNq{Q)  for  the  Lagrange  multiplier 
for  the  equality  constraint  c(9,  u)  =  0  need  to  be  chosen.  In  view  of  (1.2) 
we  suggest  A’  =  0  but  convergence  will  be  guaranteed  for  any  other  choice 
of  A’  as  well.  The  inequality  constraint  g(q)  <  0  can  be  augmented  in  a 
similar  manner  as  the  equality  constraint  e{q.u)  =  0  (see  [II\T.IK2.KKl]  for 
details).  The  hybrid  property  of  our  algorithm  is  now  evident  since  the  term 
^  I  €{q,u)  l^i  involves  the  equation  error 

I  diviqgradu)  +  /  |^-i 

and  the  term  A  |  u  —  z  |^i  represents  an  output  least  squares  criterion. 
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2  Algorithm  and  Implementation 

We  carried  out  our  computations  for  the  choices  fl  =  [0, 1]  and  Q  =  [0. 1]  x 
[0,1].  For  the  computations,  we  discretized  the  problem  using  the  finite 
element  method  [.AB];  i.e.  we  represented  q  and  u  by 

(2.1) 

i=l 

(2.2) 

J-t 

where  o"  €  and  t',’"  6  L°^(U).  We  used  the  following  types  of  basis 

functions  for  t''|"  and  d*"  in  the  two  dimensional  case. 

Type  I  The  functions  f]"  and  o’*  are  piecewise  linear  tensor  splines  (Sch): 
i.e.  m  =  (.1/  +  1 

=  Brix,)B^(x2)  for  i.j  =  0.  ■  •  • ,  .U 

and  n  =  (.V  —  1 

=  B^\xx)B^'{t2)  for  ij  =  -  1 
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where 


0  ,  elsewhere. 

Type  II  The  functions  rl'^  and  are  piecewise  linear  basis  functions 

on  triangular  elements  with  nodes  at  ,  t,j  =  0,  •  •  • ,  M  and  nodes 

at  i,j  =  ],•••, A'  —  1.  respectively,  where  m  =  {M  +  1)^  and 

n  =  (A^-  1)=. 

Type  III  The  functions  ci"  are  the  same  as  in  Type  II  and  the  functions 
V'i”  are  piecewise  constant:  i.e., 

for  i.j  =  1,  •  •  • , 

where  \{a.h)  is  the  indicator  function  of  the  interval  (a.b)  and  m  =  In 
the  one  dimensional  case  we  used  linear  splines  [Sch]  for  the  api)io.\imation 
of  both  9”‘  and  i/'’. 

As  a  regularization  term  we  took 

(2.3)  A'(?)  =  £  £  (|,„|=  +  |,„|") 

with  the  obvious  modification  in  the  one-dimensional  caise.  Here  we  did  not 
follow  the  theory  developed  in  [IK2]  which  requires  that  for  n  =2  or  3  the 
second  order  terms  be  included  and  that  N  be  chosen  as  in  (1.3). 
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With  these  specifications  made,  the  description  of  >  0,  is  given 

(FJ*’")  minimize  =  |u"  Fu" 

subject  to 

iT’^ir,u-)  ;=  (g”*  v«",  V<^?)  -  (/,<?^r )  =0  1  <  *  <  n 

qj  >  a  for  1  <  j  <  m 

where  g”  =  col(g” ,  •  •  • ,  g™J  €  is  the  coordinate  vector  of  g"*  (s’)  €  L°‘  {Ct) 
in  (2.1)  and  w"  =  col(u".  •••,«")  €  3?"  is  the  coordinate  vector  of  t/"(a’)  £ 
Ho{(l)  in  (2.2).  The  coordinate  vector  in  3?"  of  e"*’”  is  given  by  = 
fj~i^rn.n  what  follows  wc  will  usc  the  symbol  g*”  for  both  the  element 
in  L'^iQ)  and  its  vector  coordinate  in  TJ”*  and  similarly  for  i/".  Moreover 
the  bar  is  deleted  from  the  notation  of  c"*'".  The  matrices  //  €  3?"*”  and 
ir  £  and  the  vector  6  £  3?"  are  given  by 

Hi.j  =  = 

t"  =  (Vor.V2)  ?  =  l.---.n. 

We  did  not  implement  the  norm  constraint  \q\]f2  <  7  in  the  two  dimensional 
case,  see  also  section  3  in  this  respect.  If  one  were  to  implement  this  con- 
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strain! ,  one  efficient  way  is  to  choose  ir  as  the  symmetric  positive  definite 
matrix  H'  on  3?”’  given  by 

W  =  -\-W  +  Q 

with  Qj  =  (V’r,V’f)L2,  =  !,•••, m. 

The  augmented  Lagrangian  function  associated  with  (PJ*’”)  with  only  the 
equalit}’  constraint  augmented  has  the  form 

(2.4) 

The  Lagrange  multiplier  AJ  €  3i”  is  updated  by 

(2.5)  =A^  +  c*//->e-”(9r,<) 

where  Aj(j  )  =  H"-) (AJ ),o"(j)  €  Nq  approximates  the  Lagrange  multiplier 
A’  6  Hq  for  the  equality  constraint  €{q.v)  =  0  in  {P'^).  We  recall  that 

(2.5)  can  be  considered  as  a  steepest  ascend  algorithm  for  the  dual  problem 
associated  with  (PJ*'"). 

Next  we  describe  two  algorithms  to  solve  (PJ*’"). 

Algorithm  1 

(1)  Choose  Ai  =  0  and  {c*}  monotonically  increasing  with  Cj  sufficiently 
large. 
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(2)  Set  k  =  \. 

(3)  Determine  (9™, uj)  from 

minimize  AJ)  subject  to  9”*  >  a. 

(4)  If  convergence  is  achieved,  then  stop.  Otherwise,  put 

Return  to  (3)  with  k  =  k  +  1. 


The  second  algorithm  takes  advantage  of  the  quadratic  properly  of  Z’"  "; 
i.e.,  Z”'"  is  quadratic  in  q’^i  resp.  u")  if  u"  is  fixed  (resp.  if  q”'  is  fixed). 

AlgorithriT  2 

(1)  Same  as  in  Algorithm  1. 

(2)  Set  A-  =  1  and  Uq  = 

(3)  Determine  q^'  from 


(/",„)  minimize  +  Af  c’"  "(9"''.ur.,)  over  9" 


subject  to  9'"  >  a. 

(4)  Determine  wj  from 

minimize  +  +Af  €"•’"(9^,11") 

+  ^e”*'"(9?’,t/”)^//-*c’"-”(9^,w")  over  t/" 

(5)  is  identical  to  (4)  in  Algorithm  1. 
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We  point  out  that  the  choice  Uq  =  H  guarantees  that  the  solution 
g”  of  (Pe,u)  coincides  with  the  solution  of  the  equation  error  formulation  in 
H-\ 

Based  on  our  theoretical  work  [IKl.  IK2].  the  choice  of  Cj  must  be  made 
such  that  the  Lagrangian  function  of  (P^)  is  uniformly  convex  at  the  solution 
Thus  Cl  depends  on  2,/,7, /9,q  and  embedding  constants.  In  our 
calculations,  we  chose  Ci  heuristically.  The  monotonically  increasing  sequence 
{citlfcij  could  be  determined  according  to  known  heuristics  for  augmented 
Lagrangian  method,  see  [PT],  for  example.  We  took  c*  to  be  a  fixed  positive 
constant  in  most  of  our  calculations.  The  gradient  of  with  respect  to 
g”  is  given  by 

^H'g’”  +  [Xl  +  c,.p-’(^>(t/")g"’  -  /")] 

and  the  gradient  of  I”  "  with  respect  to  u"  is  given  by 

//u"  +  +  P(g'")  \Xl  +  -  .r )] 

where  e'"'”(g"'. «")  =  4>(!/")g"'  —  /"  =  H{q^)ii^  —  /"  and  the  matrices 
$(u”)  €  7^”*”’  and  H{q”')  €  are  defined  by 

^(u”)t./  =  (0r  V  V^*>»  Ar=l, •••,«;  /=!,•■ -,>77, • 

With  our  choice  of  basis  functions  for  and  the  matrices  4>(i/”)  and 
H{q"')  are  sparse.  In  the  two  dimensional  case,  the  gradient  calculations  can 
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be  performed  by  the  order  of  +  ^P  operations.  The  calculation  of 
for  given  e  €  requires  the  order  of  operations. 

Remarks  (1)  In  Algorithm  2,  Steps  (3)  and  (4)  are  a  quadratic  minimization 
in  q”'  and  u",  respectively. 

(2)  One  can  successively  use  the  Steps  (3)  and  (4)  in  Algorithm  2  to  obtain 

the  solution  to  the  minimization  of  u”, AJ)  in  Step  (3)  of 

Algorithm  1.  It  is  recommended  that  in  the  initial  stages  (i.e.,  k  =  1,2)  of 
Algorithm  2,  Steps  (3)  and  (4)  are  rei>eated  a  few  times  before  updating  the 
Lagrange  multiplier  in  Step  (5). 

(3)  VVe  used  the  conjugate  gradient  method  [AB]  to  perform  Step  (3)  in 
Algorithm  1  and  Steps  (3)  and  (4)  in  Algorithm  2  (we  used  the  ZXCGR- 
routine  in  the  IMSL  library  in  our  calculations). 

(4)  The  matrix  operation  //“*  (which  approximates  (—A)”’)  plays  the  role 
of  (pre-)  conditioning.  In  fact,  (Pout)  can  be  written  as 

minimize 


in"^(7/  +c,//(9r)//'’//(9r)K 


(2.6) 


+w"^(6"  +  P{q^)i^k  ~  +  constant  . 

Since //(9^)  is  the  approximation  to  —  V(9r(Vw))  on  V",  H  and  H(q'k)H~^H(q^) 
possess  the  same  order  of  spectral  condition  number.  Without  the  opera¬ 
tion  the  equation  error  term  in  {Pout)  would  exhibit  the  square  of  the 
condition  number  of  H.  For  the  case  H  =  [0, 1]  x  [0, 1]  (in  general, in  the 
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multi-dimensional  case)  the  operation  //"*e  for  given  e  6  3?"  is  the  most 
costly  calculation  in  our  algorithms.  Hence  one  may  use  a  smoothing  matrix 
V  based  on  the  SSOR  (symmetric  successive  overrelaxation)  method  [AB] 
or  on  a  multi-grid  method  (HT)  in  place  of  so  that  less  calculation  is 
required.  In  this  case  the  condition  number  of  H{q^)VH(q]^)  is  improved 
(see  [AB]). 

(5)  It  is  recommended  to  use  a  pre-conditioned  conjugate  gradient  algo¬ 
rithm  for  the  minimization  of  (2.6)  since  the  condition  number  of  the  matrix 
//(</r)//-’//(9r)  is  still  of  the  order  of  N^.  The  preconditioned  conjugate 
gradient  method  can  be  formulated  as  follows:  setting 

the  7-th  step  consists  of 


Q,  = 

gJhi/djQd, 

“j+1  = 

Uj  -f  Ojdj 

9}+i  — 

9,  +  0}Qdj 

hj+\  = 

9j+ihj^il9jfi 

— A7+1  +  0jdj 

where  uj  is  the  j-th  iterate  for  the  minimizer,  gj  is  the  corresponding  gradient 
vector,  and  dj  is  the  conjugate  direction.  Initially  one  chooses  uq  and  puts 
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gc  —  Huft+b" f*),  fco  =  *n<i  —  ~fto-  Again  one 

could  replace  the  smoothing  operation  by  one  of  the  alternates  described 

above. 


15 


3  Theoretical  Considerations 

In  this  section  we  will  discuss  some  of  the  theoretical  issues  with  regards 
to  estimation  problems  in  elliptic  PDEs  and  their  solution  by  means  of  the 
augmented  Lagrangian  technique. 

Convergence  of  algorithm.  The  convergence  of  the  augmented  Lagrangian 
method  that  we  described  in  the  previous  sections  has  been  studied  in 
Under  appropriate  conditions  on  the  problem  data  (2, /,  q, 

X  appearing  in  (P^),  convergence  of  the  sequence  A*) 

to  A*)  in  H^{Q)  X  Hq{Q)  x  Hl(Q)  was  shown.  Here  {qk-  Vk)  are  mini- 

mizers  (in  a  neighborhood  of  of  the  functionals  Le^(q,  u.  A*')  subject 

to  l{q)  <  0  and  ff{q)  <  0  as  specified  in  (1.4),  A^  is  determined  by  (1.5) 
and  (q‘^.tt‘^)  is  a  solution  of  (P^)  with  associated  Lagrange  multiplier  A*. 
A  critical  step  in  the  proof  of  the  convergence  results  involves  showing  the 
positivity  of  the  Hessian  of  the  Lagrangian  associated  with  (P'^),  evaluated 
at  (9^,11^,  A’).  It  is  our  future  interest  to  study  convergence  properties  of  the 
solutions  of  the  discretized  problem  (P^’")  as  {m.n)  —*  oc. 

Semi-norm  regularization.  The  use  of  the  regularization  term  3X{q)  in 
(P^)  is  common  in  solving  ill-posed  inverse  problems.  In  [CK1,CK2,KS]  for 
example,  the  stability  of  the  solutions  to  the  regularized  output  least  squares 
method  with  respect  to  perturbations  in  the  problem  data  in  (P^)  is  studied. 
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We  point  out  the  fact  that  the  regularization  term  that  we  propose  licre  is 
only  a  semi  norm  on  More  specifically,  it  does  not  contain  the  |^|i,j-term 
so  that  onij'  the  variation  of  q  is  regularized.  If  the  Zr2-norm  is  included  in 
the  regularization  term,  we  normally  obtain  under-estimated  solutions. 
Identifiability  and  Singular  set.  It  is  easy  to  argue  the  existence  of  a 
solution  of  (P^)  for  ^  >  0.  Convergence  of  (q‘^,u^)  to  a  minimum 

norm  solution  of  (P°)  (that  is,  a  solution  of  (P®)  that  minimizes  A’{<y))  as 
/3  — +  O'*'  was  studied  in  [II\2,EKN].  Uniqueness  of  the  solutions  of  (P°)  and 
(P'^),  13  >  0  cannot  be  guaranteed  in  general.  We  mention  here  a  related 
result  concerning  the  injectivity  of  the  mapping  q  —*  u{q)  and  refer  to  [K] 
for  further  discussion  of  these  matters.  Injectivity  of  q  u{q)  is  a  necessary 
but  not  a  sufiicient  condition  to  obtain  the  continuity  of  the  inverse  of  the 
map  q  it{q).  Let  q’  and  q  satisfy  the  constraints  in  (P’^)  where  q'  is  the 
“true’'  reference  coefficient  and  ^  is  a  perturbed  one.  Let  us  denote  by  u{(f} 
and  w(g),  the  corresponding  solutions.  Then  we  obtain  the  estimate  [K]. 

(3.1)  A’|u(9*)-u(g)|„.2,i(n, . 

where  the  constant  K  >  0  depends  on 

k'lu  uf-klu  '  P.  p  >  n  and  |u(9*)lx.oo. 

Let  us  define  the  singular  set 

(3.2)  5  =  {i  €  fi,  V»(?-)(x)  =  0}  . 
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Clearly,  if  meas(5)  =  0  and  u(g‘)  =  v(q),  then  q’  =  q  a.e.,  so  that  tlie 
injectivity  of  the  map  q  —*  uiq)  holds  at  q’. 

Choice  of  the  Output  Least  Squares  criterion.  In  our  calculations,  we 
use  the  /fo-norm  as  an  output  least  squares  criterion.  Such  a  choice  is  based 
on  the  facts  that  (i)  the  //o-topology  for  the  output  least  squares  criterion 
and  the  //“'-topology  for  the  equation  error  term  are  natural  from  the  point 
of  view  of  the  second  order  sufficient  optimality  condition  for  (P^)  (see  [1K2] 
for  details),  (ii)  the  choice  of  these  topologies  also  leads  to  a  method  that 
requires  the  same  amounts  of  numerical  differentiations  in  both  the  equation 
error  and  the  output  least  squares  term,  and  (iii)  it  enhances  the  sensitivity  of 
the  map  q  —*  u{q).  In  terms  of  computational  efforts  there  are  no  differences 
between  the  use  of  the  Ij'^orm  criterion  and  that  of  the  HQ-norm  criterion. 
With  a  view  on  balancing  the  output  least  squares  and  the  equation  error 
terms,  one  may  take  the  //“^-topology  for  the  equation  error  term  if  the 
norm  is  used  for  the  output  least  squares  criterion,  so  that  the  augmented 
Lagrangian  function  (1.4)  is  replaced  by 


(3.3) 


Lc^{q.v,X'’)  =  i|u  -  -I- /3.V(9) 

■\-{X,e{q,u))ii  +  ^le(g.u)lJ^j 


where  e{q,u)  =  (— A)“'(div(9  grad  u)  +  /). 
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A  general  class  of  inverse  problems  in  elliptic  partial  differential 
equations.  In  this  paper  we  restricted  our  attention  to  the  Dirichlet  bound¬ 
ary  value  problem  (1.1).  Without  difficulties  one  can  modify  our  formulation 
to  treat  other  tj'pes  of  boundary  conditions  as  well.  For  example,  consider 
the  problem  of  determining  q(x)  in 

-dW(q  grad  u)  +  p(i)u  =  / 

with  the  boundary  condition  ^  on  F  where  denotes  the  outward 

normal  derivative.  p{x)  €  L°^  with  p{t)  >  tc  >  0  and  g  €  In  this 

case,  we  define  the  equation  error  term  by 

€{q,  u)  =  (-  A  -f  l)*‘(-div{9  grad  u)  +  pu  -  /) 

where  A  is  the  Laplacian  with  dom{A)  =  {u  6  and  =  0}  and 

is  replaced  by  //’  in  While  we  are  concerned  with  estimating 

the  diffusion  coefficient  in  thi>  pai>er.  our  formulation  enables  us  to  consider 
the  estimation  problem  for  other  coefficients  in  elliptic  partial  differential 
equations  (e.g.,  p(x)  €  in  the  above  example). 

Inequality  constraints.  As  pointed  out  in  the  Introduction,  the  inequality 
constraints  (i.e.,  the  pointwise  constraint  /  and  the  norm  constraint  g)  are 
required  for  existence  of  solutions  to  (P^).  In  our  calculations,  we  ignored 
the  pointwise  lower  bound  without  harm.  Besides,  in  our  formulation  the 
positivity  of  q{T)  is  not  a  hard  constraint  since  we  do  not  need  to  solve  the 
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equation  (1.1)  for  u  in  each  iteration.  The  relevance  of  the  norm  constraint 
was  investigated  for  the  one  dimensional  c«ise  in  [KKl].  It  was  observed 
that  its  use  improved  the  results  if  the  interior  of  the  singular  set  5  =  € 

(0,  l],u,(9*)  =  0}  is  not  empty  and  it  did  not  change  them  significantly  if 
meas  5  =  0,  pro\-ided,  of  course,  that  7  is  chosen  such  that  l9*l//j(o.])  £  7- 
We  did  not  use  the  norm  constraint  in  our  implementation  for  the  case  Q  = 
(0, 1]  X  [0, 1]  since  we  do  not  expect  that  any  new  phenomena  arise.  Moicowr 
the  regularization  term  also  provides  a  semi-norm  bound  for  the  solutions. 
Miscellaneous  comments. 

An  implicit  regularization  is  often  achieved  by  choosing  coarser  basis  el¬ 
ements  for  representing  q{x)  than  for  w(x).  For  example,  in  some  of  onr 
calculations,  we  took  .V  =  2.1/  (see  (KKl,  KK2.  KK3]). 
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4  Numerical  Results 


Numerous  numerical  experiments  were  carried  out  for  the  two  algorithms 
that  we  explained  in  Section  2.  The  tests  were  conducted  in  the  following 
w’ay.  For  all  examples  except  for  Example  7  we  chose  the  ‘*true“  parameter  q’ 
as  well  as  the  solution  u’  =  u{q“)  and  calculated  /  from  (1.1).  For  Example 
7  we  chose  q‘  and  /  and  calculated  «*  =  u(q‘)  numerically.  The  data  r  were 
determined  from  u{q‘)  either  by  putting  z  =  v{q“)  (“distributed  data")  or 
by  evaluating  u(q')  at  points  and  determining  ?  as  a  cubic  Hermite 

spline  interpolation  of  z,  =  u(q’.T,).  Noisy  data  were  obtained  by  specifying 

z,  =  t/(?‘,x,)  +  6^,, 

with  uniformly  distributed  random  number  in  [-1,1]^  and  choosing  r  as  the 
cubic  Hermite  spline  interpolation  of  these  z,.  The  numerical  examples  to 
be  discussed  below  are  divided  into  two  groups  an  are  based  on  iwo  different 
packages.  The  first  six  examples  demonstrate  the  performance  of  the  algo¬ 
rithm  while  a  specific  aspect,  as  for  instance  the  behavior  near  the  singular 
set,  is  under  investigation  and  the  number  of  unknowns  is  no  larger  than  SI. 
In  Examples  7  and  8  we  show  that  the  algorithm  that  we  propose  also  works 
very  well  for  an  extremely  fine  resolution  (C4  x  64)  of  the  unknown  coeffi¬ 
cient.  For  the  calculations  of  Examples  1-6  we  took  Type  I  basis  functions 
with  the  grid  for  the  state  approximation  twice  as  fine  as  the  grid  for  the 
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coefficient  approximation  (N  =  2Af).  This  choice  introduces  some  inherent 
regularization.  The  (— A)“*  operation  is  resolved  bj'  Cholesky  factorization 
and  the  optimization  problems  are  solved  with  the  IMSL  version  ZXCGR 
of  the  conjugate  gradient  algorithm.  The  discretization  for  Examples  7  and 
8  used  Type  III  basis  functions  with  N  =  M.  The  (  — A)“*  operation  is  real¬ 
ized  by  a  multigrid  step  involving  one  V-cycle  with  finite  element  injection, 
w’hile  the  optimization  problems  are  solved  by  a  preconditioned  conjugate 
gradient  algorithm  as  explained  in  (5)  of  Section  2.  Moreover,  Algorithm  2 
was  slightly  modified  by  iterating  several  times  between  steps  (3)  and  (4) 
before  going  to  step  (5).  The  results  from  a  coarser  mesh  are  then  used  as 
start-up  values  for  the  next  finer  mesh.  The  calculations  for  Examples  1-0 
were  carried  out  on  an  IBM- AT  and  those  for  Examples  7-8  on  an  IBM  3081. 
both  in  double  precision.  For  further  specification  of  the  numerical  aspects 
and  for  additional  test  examples  we  refer  to  [KK  1-3]. 

W’e  summarize  some  general  observations. 

(1)  In  all  examples  the  augmented  Lagrangian  algorithm  in  the  form  of 
Algorithm  1  as  well  as  Algorithm  2  performed  well. 

(2)  In  all  examples  where  complete  data  were  available  (i.e.  z  =  u{q')). 
where  q’  was  a  smooth  function  and  where  the  singular  set  5  w£is  small,  the 
first  or  second  iteration  gave  a  good  approximation  to  both  q'  and  u{q‘). 
If  these  requirements  were  not  satisfied,  and  m  was  sufficiently  large,  then 


higher  iterations  gave  an  improvement  over  the  first  one.  Recall  that  for 
Algorithm  2,  the  determination  of  the  first  iteration  coincides  with  the 
solution  of  the  equation  error  algorithm. 

(3)  Algorithm  2  is  significantly  faster  than  Algorithm  1. 

(4)  For  the  estimation  of  q’  the  form  of  the  singular  set  S  is  of  great 
importance.  If  the  dimension  of  S  is  smaller  than  the  dimension  of  Q,  (and 
/3  =  0  OT  0  is  very  small),  then  the  maximum  error  of  the  approximation  to 
q‘  generally  occurs  within  a  small  neighborhood  of  5.  If  S  contains  an  open 
set,  then  q'  cannot  be  estimated  there. 

(5)  We  could  observe  conversion  and  rates  of  convergence  of  the  solutions 
to  (/’J''")  as  m  — >  oc,  n  — *  oc,  see  (KKI,  KK2]. 

(6)  The  use  of  a  regularization  term,  especially  if  the  same  grid  is  cho¬ 
sen  for  both  the  coefficient  and  the  state  space  approximation,  or  to  avoid 
oscillations  in  the  neighborhood  of  a  singular  set,  is  useful.  We  are  currently 
working  on  an  adaptive  algorithm  to  specify  the  level  of  regularization  that 
is  required  for  a  specific  problem. 

(7)  The  augmented  Lagrangian  algorithm  that  we  propose  also  works 
very  well  when  applied  to  problems  with  discontinuous  unknown  coefficients. 
Comparing  the  results  where  a  gridpoint  of  the  discretization  for  the  coeffi¬ 
cient  coincides  with  a  discontinuity  of  g*,  to  those  where  discontinuities  and 
grid  points  do  not  coincide,  the  latter  were  clearly  better. 
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Most  of  the  plots  that  we  present  below  are  endowed  with  the  following 
information: 

Example  X,  0  =  X,  M  =  X,  k  =  A',  L*  =  x,  additional  comments. 
Here  stands  for  the  i^-error  Unless  otherwise  specified, 

the  data  were  taken  to  be  error  free  distributed  observations,  2  = 
Examples  3-6  were  calculated  with  Algorithm  1,  the  remaining  examples 
with  Algorithm  2.  Since  the  numerical  results  always  produced  a  very  good 
approximation  to  u{q‘)  we  only  show  the  numerical  results  for  the  approxi¬ 
mation  of  the  coefficient  q"  (except  in  Example  6).  In  the  one  dimensional 
examples  =  (0, 1)  and  in  the  two  dimensional  examples  ft  =  (0, 1)  x  (0, 1), 
The  start-up  value  for  q”^  for  the  finite  dimensional  approximating  problems 
(pm.ri)  chosen  to  be  identically  1  in  Examples  1-6  and  1.5  in  Examples 
7  and  8.  For  all  calculations  c*  weis  taken  to  be  equal  to  1. 

Example  1.  Here  we  take 

2(.T,j/)  -  u{q'){x,y)  =  uix)u'{y), 

q'  =  2-1- sin(j-V) 

where 

I  -9x*  -I-  6x  ,  X  €  [0, 

I  -9x2-|-12x-|.3  ,x6[f,l], 

with  /  calculated  from  (1.1).  In  Plot  1  we  give  the  graph  for  2  and  q*  as  well 
as  the  numerical  result  for  q  after  8  iterations  without  and  with  the  use  of  a 
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regularization  term.  Here  the  singular  set  has  non-empty  interior.  Observe 
that  regardless  of  the  use  of  a  regularization  term,  g*  is  identified  well  over 
the  complement  of  the  singular  set. 

Example  2.  Here  we  took 

z{x,y)  =  u(g’)(a:,y)  =  (y  -  2y®)sin27rxsin2ry 
g'  =  2  -f  sin(x^y). 

and  /  is  calculated  from  (1.1). 

The  singular  set  for  this  example  consists  of  isolated  points  and  the  line 
characterized  by  y  =  0  and  y  =  \.  Plot  2  gives  a  graph  of  |  V  x|  and  Plot 
3  shows  the  results  after  the  first  and  eighth  iteration.  We  point  out  the 
improvement  in  the  neighborhood  of  the  singular  set  characterized  by  y  =  5. 
Example  3.  Here  we  took 

t/(g*)  =  sin  27rx  sin  'lity 

g*  -  1 -f6xV(l -y). 

and  calculated  /  from  (1.1).  Then  tests  were  carried  out  assuming  that  suc¬ 
cessively  more  data  points  become  available.  The  results  after  one  iteration 
of  Algorithm  1  are  shown  in  Plot  4.  Further  iterations  so  not  change  tlic 
results  significantly  for  this  example.  This  is  probably  due  to  the  fact  that 
the  discretization  of  g  is  rather  coarse. 
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Example  4.  This  is  an  example  with  a  discontinuous  coefficient; 

r  =  u(9*)  =  sinffx 

[  1  /or  X  6  [0,  i] 

9’  =  ■{  2-x  /orx€(|,f) 

[  — 9x^  +  yx  —  3  for  X  6  [f,  1]. 

The  results  for  M  —  21,  M  =  22  and  M  =  42  are  given  in  Plot  5.  For  M  =  21 
and  M  =  42,  the  discontinuities  are  also  grid  points  of  the  discretization  for  q. 
We  point  out  that  these  estimates  are  obtained  without  a-priori  assumptions 
on  the  location  of  the  discontinuity  and  the  size  of  the  jump.  The  results  can 
be  improved  by  using  regularization,  see  Plot  6. 

Example  5.  Here  q‘  is  chosen  as  in  Example  4.  We  compare  the  nu¬ 
merical  results  for  two  different  observations  Zi  =  sinrx  and  *2  =  sin2;rx. 
Observe  that  the  singular  set  for  Zt  is  5(2i)  =  ^  and  similarly  5(22)  =  f  }• 

The  last  graph  in  Plot  7  shows  the  numerical  result,  when  .Algorithm  1  was 
modified  in  the  obvious  way  to  account  for  two  observations  simultaneously. 
Example  6.  For  this  example  we  chose 

u{q")  =sin7rx. 
q‘  =  1  +  X. 

We  put  z{xi)  =  u(9')(xi)  with  x,  =  for  i  =  0,'--,7,9, •••,42,  and 
2(xg)  =  u(q’){Xi)  +  J,  thus  producing  an  outlier  at  Then  w’e  passed 
a  cubic  Hermite  spline  interpolate  through  Zj,  see  the  thicker  lines  in  Plot 
8.  To  estimate  the  unknown  parameter  9*  we  used  Algorithm  1  as  well 
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as  Algorithm  1  with  the  [u  —  z\if\  output  criterion  replaced  by  an  |u  — 
criterion.  The  second  curve  in  the  graphs  of  Plot  8  give  the  numerical 
result  for  u".  Differently  from  the  “’-criterion,  the  outlier  is  essentially 
ignored  by  the  numerical  solution  u**,  when  the  -criterion  is  used.  In 
this  Ctise  the  approximation  q”'  to  9*  is  still  qualitatively  correct  while  the 
//’  criterion  without  regularization  produces  a  useless  estimate  for  q’ .  For 
further  discussion  of  this  example  see  {KK3]. 

Example  7.  Here  we  took 

/  =  sin  2kx  sin  2jry, 

f  2  /or{x,y)e(.3,.3)x{.3,.6) 

^  ~  \  1  elsewhere 

and  data  were  assumed  to  be  available  on  a  20  x  20  grid  for  the  first  result 
and  on  a  40  X  40  grid  for  the  second  result.  Observe  that  adding  more  data 
points  results  in  a  sharper  resolution  of  the  discontinuous  coefficient  q‘. 
Example  8.  For  this  example,  with 

u(9* )  =  sin  2xx  sin  2xy, 
q’  =l+6x^y{l-y) 

consecutively  more  noise  was  added  to  u(9")(a‘,)  at  the  observation  points  a;, 
which  were  assumed  at  a  uniform  10  x  10  grid.  Observe  that  we  are  referring 
to  absolute,  rather  than  relative  noise  here. 
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Example  1  Example  1 

function  z  function  q’. 

Plot  1 

Example  1  Example  1 

=  0,  A/  =  6,  it  =  8,  I*  =  0.12  ^  =  0.  A/  =  6,  it  =  8,  = 

Plot  2  Example  2 
^  function  1  V  si 

PI  .  _  Example  2  Example  2 

*  ^  =  0,  A/  =  9,  it  =  1,  £2  =  0.034  ^  =  0,  A/  =  9,  it  =  8,  = 

Example  3  Example  3 

Plot  4  =  0,  A/  =  5,  it  =  1  ^  =  0,  A/  =  5,  it  =  1 

interpolated  data  at  {|,^}  interpolated  data  at 
Example  3 

;?  =  0,  A/  =  5,  it  =  1 
interpolated  data  at 

Example  4  Example  4 

=  0,  A/  =  21.  k  -2,  =  0.101  13  =  0,  M  =  22,  k  =  ,  P 

Plots 

Example  4  Example  4 

/?  =  0,  A/  =  41.  k  =  2,  £2  =  0.055  ^  =  0,  M  =  42,  k  =  , 

Pint  fi  Example  4 

^  a  =  10-®,  M  =  41.  k  =  3,  £2  =  0.0445 

Example  5  Example  5 

a  =  0,  M  =  41,  A'  =  3,  £2  =  0.0525  /?  =  0,  A/  =  41,  A*  =  3. 

Example  5 

^  =  0,  A/  =  41,  A-  =  3,  £2  =  0.0447, 
two  distributed  observations 

PI  .  j.  Example  6  Example  6 

/?  =  0,  A/  =  6,  A'  =  3,  W’’ •^-criterion  =  0.1.  M  =  6,  k  =  3, 

Example  7 

a  =  10“^,  A/  =  64,  k  =  20  X  20  point  measurements 

Plot  9 

Example  7 

a  =  10“^,  M  =  64,  A:  s=  40  x  40  point  measurements 


0.0036 

=  0.0065 

=  0.064 

=  0.060 

£2  =  0.017 

IP  -criterion 
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Plot  10 

/3 


Example  8 

0  =  10“*,  M  =  64,  it  =  10  X  10  point  data  with  0%  noise 
a=  5  X  10“*,  A/  =  64,  it  =  10  X  10  point  data  with  1%  noise 
=  5  X  10“*,  A/  =  64,  it  =  10  X  10  point  data  with  Z%  noise 
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